Tracing ultrafast interatomic electronic decay processes in real time and space 
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Tremendous advances in laser pump-probe techniques open the door for the observation in real 
time of ultrafast electronic processes. Particularly attractive is the visualization of interatomic 
processes where one can follow the process of energy transfer from one atom to another. The inter- 
atomic Coulombic decay (ICD) provides such a process which is abundant in nature. A multielectron 
wavepacket propagation method enables now to trace fully ab initio the electron dynamics of the 
process in real time and in space taking into account all electrons of the system and their corre- 
lations. The evolution of the electronic cloud throughout the ICD process in the rare gas cluster 
NeAr following Ne2s ionization is computed and analyzed. The process takes place on a femtosecond 
timescale, and a surprisingly strong response is found at a much shorter attosecond timescale. 

PACS numbers: 34.30.+h, 31.70.Hq, 82.33.Fg 



In recent years tremendous developments of laser 
pump-probe experimental techniques made possible a di- 
rect observation in real time of different kinds of ultra- 
fast processes with sub- femtosecond resolution [H, [H) H 0] ■ 
This opened the door for the investigation of electronic 
processes in the attosecond/femtosecond timescale that 
take place before the nuclear dynamics comes into play. 
Processes like the rearrangement of the electronic system 
following an excitation of an inner-shell electron can now 
be traced in time and space and analyzed. These capabil- 
ities promise a revolution in our microscopic knowledge 
and understanding of matter. 

The excitation of an electron from an atomic shell other 
than the outermost valence orbital creates a transient 
hole state. This hole is not stable and the system tends 
to minimize its energy by filling the vacancy with an elec- 
tron from an outer shell, the excess binding energy be- 
ing carried away either by an extreme UV or X-ray flu- 
orescence photon, or transferred via Coulomb forces to 
another electron, which subsequently escapes from the 
atomic binding. When possible energetically, the latter 
non-radiative mechanism of de-excitation is extremely ef- 
ficient in comparison to the competing photoncmission. 
Indeed, the characteristic times are typically less than 
100 femtoseconds (1 fs = 10~ 15 s) to even well less 
than a femtosecond entering the attosecond regime (1 as 
= 10~ 18 s), compared to radiative decay lifetimes which 
except for core levels of heavy elements, belong to the 
nanosecond range (1 ns = 10~ 9 s). 

The non-radiative decay processes can be divided into 
two major categories depending on whether the electrons 
involved in the process belong to the same or to different 
subunits of the system. The former are referred to as 
intra-atomic/molecular decay processes and the latter as 
inter- atomic/molecular decay processes. An example of 
intra-atomic decay is the well known Auger effect, follow- 
ing a core ionization of atoms, which besides its funda- 
mental importance has a wide range of applications. De- 



pending on the energy of the core hole and the strength 
of the electronic coupling, the Auger lifetimes range from 
few femtoseconds to few hundreds attoseconds (in the so- 
called super-Coster-Kronig transitions). That is why the 
observation in real time of the electron dynamics of such 
processes became possible only recently. Few years ago, 
Drescher et al. [l[ using a sub-femtosecond X-ray pulse 
for excitation and a few-cycle light pulse for probing the 
emission traced the electron rearrangement during the 
M4 j 5NiN2 i 3 Auger decay in Kr, with lifetime of the M- 
shell vacancy of about 8 fs, giving birth to time-resolved 
atomic inner-shell spectroscopy. 

Contrary to core ionization, the ionization of inner- 
valence electrons usually produces ions in excited states 
lying energetically below the second ionization threshold 
thus making the slow radiative decay the only possible 
de-excitation mechanism as long as the resulting ions are 
isolated. It has been shown recently [5[ that the situa- 
tion is fundamentally different if the ions are embedded 
in an environment or have neighbors like in a cluster. 
Then, the possibility to distribute the positive charges 
between cluster subunits substantially lowers the double 
ionization threshold of the cluster compared to that of the 
isolated subunit giving rise to an inter-atomic/molecular 
decay mechanism, where the excess energy of the ion 
with the inner-valence hole is utilized to ionize a neu- 
tral neighbor. The process is ultrafast, in the femtosec- 
ond timescale, and the emitted electron has a kinetic en- 
ergy of a few electronvolts. This process, called inter- 
atomic (or intermolecular) Coulombic decay (ICD), was 
predicted theoretically and shown to follow a general de- 
cay mechanism, taking place both in hydrogen bonded 
and in van der Waals clusters @, @, 0|- Very recently, 
the theoretical findings have been confirmed in a series 
of spectacular experiments carried out by several groups 
[1, 0, 0, U > 13] • Other experimental groups are engaged 
in the preparation of time-resolved ICD experiments us- 
ing novel attosecond pulse techniques lj, |l4{ . 
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The goal of the present work is to provide for the hrst 
time a theoretical description of an interatomic decay 
process in real time and space. Such a description re- 
quires accurate ab initio calculations of the time evolu- 
tion of the electronic cloud including explicitly the corre- 
lations among all the electrons, i.e., it requires multielec- 
tron wavepacket dynamics. Here, we present such a com- 
putational method and apply it for tracing in time and 
space the ICD process in NeAr following the 2s ioniza- 
tion of Ne. It should be noted that the method is equally 
suitable for calculating in real time all kinds of ultrafast 
electron relaxation processes following an ionization of a 
system. The technical details of the method are given 
elsewhere [Hj], where it is used to compute electron dy- 
namics of non-decaying states. Here, only the theoret- 
ical foundations of the method will be sketched as well 
as some subtleties concerning its application to decaying 
states. We stress that treating the dynamics of a decay 
process represents a much higher degree of complexity 

The starting point of our investigation is a neutral sys- 
tem l^o)- The ionization of the system generates a non- 
stationary state The resulting hole charge then is 
traced in time and space, i.e., the time-dependent hole 
density is calculated. For convenience we assume that the 
initial state is created by the sudden removal of an elec- 
tron. The hole density of the ionized system is defined 
by 

Q(r,t) := (*o|p(»',*)|*o>-(*<|p(» ! ;t)|*<>=Po(» s )-Pi(» ! ;t) J 

where p is the local density operator, and is the gen- 
erated initial cationic state. The first term in Eq. ((TJ 
is the time-independent ground state density of the neu- 
tral system, p$, and the second one, Pi, is the density of 
the cation which is time-dependent, since is not an 
eigenstate of the cation. The quantity Q(r,t) describes 
the density of the hole at position r and time t and by 
construction is normalized at all times t. In the Heisen- 
berg picture, the time-dependent part pi(r,t) reads: 



(2) 

where |3>i(t)) = e _ ^^ fft |$i) is the propagating multi- 
electron wavepacket. 

Using the standard representation of the density oper- 
ator in a one-particle basis {<p P (r)}, often called orbitals, 
and occupation numbers {n p }, Eq. (TTJ) can be rewritten 
in the following form 



Q(f,t) = £>;(^(r)AW(*), 



(3) 



where the matrix N(t) = {N pt] (t)} with elements 
N pq (t) = 6 pqTlp 

M.N 



is referred to as the hole density matrix. The second term 
of Eq. ([4J is obtained by inserting in Eq. ([2]) a resolution 
of identity of a complete set of appropriate ionic eigen- 
states I^m)- The matrix pmn is the representation of 
the density operator within this basis. 

Diagonalization of the matrix N(t) for fixed time 
points t leads to the following expression for the hole 
density 



Q(r,t) = ^2 \(p p {r,t)\ 2 n p {t), 



(5) 



where <p p (r, t) are called natural charge orbitals, and n p (t) 
are their hole occupation numbers. The hole occupation 
number, n p (t), contains the information which part of 
the created hole charge is in the natural charge orbital 
tp p (r,t) at time t. Because of the conservation of hole 
charge, one finds that J2 P n P (t) = 1 at any time. The hole 
occupation numbers, together with the hole density, are 
central quantities in the observation and interpretation of 
the multielectron dynamics taking place after the removal 
of an electron. 

For calculating the hole density matrix, Eq. ((4J, we 
have used ab initio methods only. The description of the 
non-stationary ionic state was done by means of the for- 
malism of Green's functions, using the so-called alqebraic 
tiagrammatic construction (ADC) scheme 00, and 
the direct time propagation of the electronic wavepacket 
was performed through the short iterative Lanczos tech 
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nique (see, e.g., Ref. |18jl)- For more details see Ref. 
and references therein. An important point should be 
addressed. Since the formalism is used for tracing the 
evolution of decaying states, i.e., a second hole is created 
in the system and an electron is ejected into the con- 
tinuum, special care must be taken in constructing an 
appropriate basis set in order to have a good description 
of the continuum electron at least in some volume around 
the origin. 

We have applied the above sketched methodology to 
describe in real time the interatomic decay of the Ne2s 
vacancy in NeAr. The choice of the system is motivated 
by the availability of Ne„Ar. m clusters [Hj], an( i by the 
fact that the lifetime of the Ne2s hole (35fs, see below) 
is short compared to the nuclear motion (vibrational pe- 
riod of NeAr is 1.2 ps and the rotational period is 180 
ps), justifying the use of clamped nuclei. Furthermore, 
the presentation simplifies by choosing a heteroatomic 
system where the initially ionized atom (Ne) is well dis- 
tinguishable for its neighbor (Ar). The internuclear dis- 
tance is taken to be 3.5 A, the equilibrium geometry of 
the NeAr cluster. To simulate an experiment one would 
have to compute the process at various internuclear dis- 
tances and then to average over the weighted distribution 
of these distances. Such a procedure is too costy in view 
of the very large effort already invested in the present 
calculations. Nevertheless, such an averaging should only 
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induce some smoothing of the curves presented below and 
a small shift of the predicted decay time. 

The electron dynamics calculations on NeAr have been 
performed using a combination of atomic and distributed 
Gaussian basis sets. The atomic basis set was chosen 
to be d-aug-cc-pVDZ on Ne and aug-cc-pVDZ on Ar. 
The distributed Gaussian basis consisted of (6s,4p,3d) 
Kaufmann-Baumeister-Jungen [l9( diffuse functions cen- 
tered between the Ne and Ar, and 36 s-functions placed 
on a grid on the x, y and z axes. The positions and the 
exponents of the latter are optimized to approximate up 
to 15 A from the system the radial part of the Coulomb 
s-wave with a kinetic energy of the ICD electron. 
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FIG. 1: (color online) Time-dependent hole occupation num- 
bers of NeAr after Ne2s ionization. Note the different 
timescale on the left part of the graphic. The initially ionized 
Ne2s orbital, plotted in black, loses more than 80 % of its hole 
charge in about 35 fs. In the same time, two other holes are 
opened - one on Ne2p (red curve), and one on Ar3p (green 
curve) - and more than 90 % of an electron is created in the 
continuum (blue curve). 

In Fig. 1 we show the results for the hole occupation 
numbers (see Eq. ([5])) as a function of time for NeAr af- 
ter sudden removal of a Ne2s-electron. This implies that 
fiNe2s = 1 at i = and all other hole occupation numbers 
are zero. As time proceeds, the initial hole (black curve) 
is gradually filled up, and two other holes are opened (red 
and green curves) . This is accompanied by the creation of 
an electron ( "negative" hole) in the virtual orbital space 
(blue curve) . For transparency, the hole occupations cor- 
responding to the p x , p y and p z components of the or- 
bitals are grouped together, as well as all negative occu- 
pations contributing to the description of the continuum 
electron. The time evolution of the hole occupations re- 
flects the timescale of the ICD process. After about 35 fs, 
the initial hole on Ne2s is filled by an electron from the 
Ne2p orbitals and an electron from the Ar3p orbitals is 
ejected into the continuum, represented in our treatment 



by a vast number of virtual orbitals. The ICD lifetime 
is thus about 35 fs, in a very good agreement with the 
result obtained by elaborate ab initio calculations of the 
decay width [2(| • Several decay channels participate (see 
below) and, therefore, the shape of the Ne2s hole occu- 
pation curve is not purely single-exponential, but rather 
a linear combination of several exponential functions. In 
principle, one may think of another process - the hole on 
Ne2s is filled by an electron from the Ar3p and an elec- 
tron from the Ne2p is ejected. However, our numerical 
analysis of the propagating electronic wavepacket shows 
that the decay probability of this pathway is negligible 
being several orders of magnitude smaller than that of 
the ICD process. 

It is worth mentioning that due to the finite basis set 
used, after some time the hole occupation numbers cease 
to reflect the physical reality when dealing with decaying 
states, since the so described continuum electron actu- 
ally cannot leave the system. When the ICD electron 
reaches the spatial end of the basis set it can be reflected 
back yielding unphysical oscillations of the hole occupa- 
tion numbers. In the present study such oscillations ap- 
pear after about 50 fs. That is why, the size of the space 
covered by the basis set is of crucial importance for the 
proper description of decaying states. 

An interesting phenomenon is observed in the ultra- 
short timescale after the sudden ionization. On the left- 
hand side of Fig. 1 the first 100 as of the process are 
presented on an expanded scale. Besides the fast drop of 
the initial occupancy for about 50 as, shown to be uni- 
versal for multielectron systems and related to the filling 
of the exchange-correlation hole 2l| , one observes an ex- 
tremely strong response of the system to the creation of 
the Ne2s hole. In just 30 as, more than half of a full 
hole is opened on Ne2p and more than half of an elec- 
tron is already in the continuum. In such a short time 
the system is already "prepared" for the consecutive ICD 
process which is completed 35 fs later. The removal of 
a Ne2s-electron is seen to introduce an enormous distur- 
bance of the electronic cloud yielding an extremely fast 
hole-particle excitation. We mention that at such short 
times a local "violation" of the energy conservation is 
possible. Indeed, following the Heisenberg uncertainty 
relation, AEAt ~ %, one finds for times ~ 30 as an en- 
ergy dispersion AE ~ 220 eV! This is by far more than 
the energy involved in the whole ICD process. 

More insight into the process is gained by inspecting 
the time evolution of the hole density, Eq. |(SJ). In Fig. 
2 the charge Q(z,t) obtained by integrating Q(r,t) over 
the coordinates perpendicular to the Ne-Ar axis (z axis) 
is plotted. The displayed results support the conclusion 
that the dominant decay is to Ne + (2pJ 1 )Ar + (3pJ 1 ); see 
also Ref. [20j. Indeed, the initial hole density (left-hand 
side of the surface in Fig. 2) changes its character from 
an one- to a two-ridge surface, corresponding to the fill- 
ing of the initial Ne2s-hole predominantly by a Ne2p z - 



4 



hole density Q(z,t) 




FIG. 2: (color online) Hole density Q(z,t) plotted against 
the Ne-Ar axis 2 as a function of time after Ne2s ionization 
of NeAr. Q(z,t) is obtained by integrating the total charge 
over the coordinates perpendicular to the z axis. The initial 
hole density changes its character from an one- to a two-ridge 
surface, corresponding to the filling of the Ne2s-hole predom- 
inantly by a Ne2p z -electron. In the same time a Ar3p-hole is 
created (right-hand side of the plot), which also has mainly 
p z -character up to about 40 fs. One clearly sees that the dom- 
inant ICD-channel is Ne + (2s" 1 )Ar^Ne + (2pJ 1 )Ar + (3pJ 1 ). 

electron. The second hole opened on Ar (right-hand side 
of the surface) displays predominantly p z -character up to 
about 40 fs. The channel ending with Ar + (3p~y), seen 
as a tiny, steadily increasing ridge situated between the 
two wide humps on the right-hand side of the surface, 
becomes important after that time indicating a slower 
decay channel. 

The surface in Fig. 2 represents the total charge den- 
sity obtained by summing over both positive and negative 
occupancies, i.e., over hole and electron charges. Because 
of the mutual cancellations in some spatial regions, it is 
insightful to separate the hole and electron charges, i.e., 
to show separately the positive and negative occupancies. 
In Fig. 3 the so determined hole charge density Q(f, t) is 
shown in 3D for the three different times t = 0, 10 and 35 
fs. One clearly sees the stages of the ICD process. After 
10 fs the initial 2s hole on Ne is already a mixture of s- 
and p z -character, and at the same time the p z -hole on Ar 
is already partly formed. After 35 fs the decay via the 
faster channel of the ICD process is essentially completed 
and we see Ne2p z and Ar3p z holes, while the Ar3p xy or- 
bitals have just started to contribute to the formation of 
the Ar3p-hole. For completeness, the electron density in 
the continuum at t = 35 fs is also displayed in Fig. 3. 

Let us summarize. A fully ab initio methodology for 
calculating in real time ultrafast electron relaxation dy- 
namics processes after ionization of a system is proposed. 
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FIG. 3: (color online) The contribution of the holes (n p (t) > 
in Eq. ([5])) to the hole density Q(r,i) at times 0, 10 and 35 
fs. For completeness, at t — 35 fs the electron density in the 
continuum (n p (t) < in Eq. (J5])) is also given (light blue). 
The colors used correspond to those utilized in Fig. 1. 

The method gives the possibility to trace in time and 
space the electron dynamics throughout different non- 
radiative decays taking into account all electrons of the 
system and their correlations. We have presented the 
first study of the multielectron dynamics of the inter- 
atomic Coulombic decay in NeAr following Ne2s ioniza- 
tion. The results disclose many microscopic details of the 
decay process in space as time proceeds. To trace the 
dynamics of the decay directly in the time domain with 
attosecond resolution, a pump-probe experiment can be 
done in analogy to that described in Ref. [l|. We hope 
that the present research will stimulate future experimen- 
tal and theoretical studies. 
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* e-mail: alexander.kuleff@pci.uni-heidelberg.de 
[1] M. Drescher et al, Nature 419, 803 (2002). 
[2] H. Niikura et al, Nature 417, 917 (2002). 
[3] H. Niikura et al., Nature 421, 826 (2003). 
[4] H. Niikura, D.M. Villeneuve and P.B. Corkum, Phys. 
Rev. Lett. 94, 083003 (2005). 



■5 



[5] L.S. Cederbaum, J. Zobeley and F. Tarantelli, Phys. Rev. 

Lett. 79, 4778 (1997). 
[6] R. Santra and L.S. Cederbaum, Phys. Rev. Lett. 90, 

153401 (2003). 

[7] V. Averbukh and L.S. Cederbaum, Phys. Rev. Lett. 96, 
053401 (2006). 

[8] S. Marburger et al, Phys. Rev. Lett. 90, 203401 (2003). 

[9] T. Jahnke et al, Phys. Rev. Lett. 93, 163401 (2004). 
[10] G. Ohrwall et al, Phys. Rev. Lett. 93, 173401 (2004). 
[11] S. Barth et al, Phys. Chem. Chem. Phys. 8, 3218 (2006). 
[12] Y. Morishita et al, Phys. Rev. Lett. 96, 243402 (2006). 
[13] S. Leone and D. Neumark, private communication. 
[14] H. Kapteyn, private communication. 
[15] A.I. Kuleff, J. Breidbach and L.S. Cederbaum, J. Chem. 



Phys. 123, 044111 (2005). 
[16] L.S. Cederbaum, in Encyclopedia of Computational 
Chemistry, P.v.R. Schleyer et al. (eds), (Wiley, Chich- 
ester, 1998). 

[17] J. Schirmer, A.B. Trofimov and G. Stelter, J. Chem. 

Phys. 109, 4734 (1998). 
[18] C. Leforestier et al, J. Comp. Phys. 94, 59 (1991). 
[19] K. Kaufmann, W. Baumeister and M. Jungen, J. Phys. 

B 22, 2223 (1989). 
[20] S. Scheit et al, J. Chem. Phys. 124, 154305 (2006). 
[21] J. Breidbach and L.S. Cederbaum, Phys. Rev. Lett. 94, 

033901 (2005). 



